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We give heuristic arguments and computer results to support the hypothesis that, after appro- 
priate rescaling, the statistics of spacings between adjacent prime numbers follows the Poisson 
distribution. The scaling transformation removes the oscillations in the NNSD of primes. These 
oscillations have the very profound period of length six. We calculate the spectral rigidity A3 for 
prime numbers by two methods and we find the cross-overs in their behaviors. 
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I. INTRODUCTION 

The primes numbers often provided a toy model for 
some physical ideas in the past. For example in [1] the 
multifractal formalism was applied to prime numbers, in 
[2] the appropriately defined Lyapunov exponents for the 
distribution of primes were calculated numerically. In 
the paper [3] it was shown that the distribution of prime 
numbers displays the 1// noise, while in [4] the noise 1//^ 
was found in the difference between the prime-number 
counting 7r{x) function and Riemann's function R{x). In 
[5] and [6] random walks on primes numbers were defined. 
In [7] an attempt to construct the dynamical model for 
prime numbers was taken and computable information 
content as well as entropy information of the set of prime 
numbers were calculated. 

The prime numbers can be regarded as eigenvalues of 
some quantum hamiltonian. The problem of construction 
of a simple one-dimensional Hamiltonian whose spectrum 
coincides with the set of primes was considered in [8] , [9] , 
[10], see also review [11]. Then it is natural to investigate 
the spacings between prime numbers, i.e. in physical lan- 
guage the nearest neighbor spacing distribution (NNSD). 
Several authors have undertaken a study of this problem 
in the past, see [12], [13], [14]. Below we will treat prime 
numbers as the energy levels and we will apply methods 
used to describe statistical properties of discrete spectra. 
Let the quantum system possess the discrete spectrum 
Ei,E2, ... and let N{E) = ©(^ " ^n) (© is a unit 
step function) denote the function counting the number 
of energy levels < E. Usually spectral staircase N{E) 
can be split into the "smooth" N{E) and fluctuating (os- 
cillating) N{E) parts. For example, for a large class of 
differential operators on d dimensional bounded manifold 
O c the Weyl's law 



N{E) 



vol(l^) 
(27r)^ 



(1) 



holds, see e.g. [15, Ch.l] . 

Given the spectrum £^1 , £^2 , • • • the statistics of nor- 
malized and dimensionless ("unfolded" spectrum, see e.g. 
[16, Sect. 4. 7]) gaps between two consecutive energy lev- 
els Sn = (£n+i — En)) /d{E)^ whcrc d{E) is the mean 
distance between energy levels up to £, was extensively 
studied in the past. For general systems £n+i — E^ are 
arbitrary real numbers and histogram of the level spac- 
ings Sn is built. It is well known, that level-spacing s 
distributions of quantum systems can be grouped into 



a few universality classes connected with the symmetry 
properties of the hamiltonians: Poisson distribution (i.e. 
P{s) = e~^) for systems with underlying regular classi- 
cal dynamics, Gaussian orthogonal ensemble (GOE, also 
called the Wigner-Dyson distribution) — hamiltonians 
invariant under time reversal, Gaussian unitary ensemble 
(GUE) — not invariant under time reversal and Gaussian 
symplectic ensemble (GSE) for half-spin systems with 
time reversal symmetry. There are many reviews on these 
topics, we cite here [17], [16], [18]. 

There is some confusion regarding the proper statistics 
of the gaps between consecutive primes: in [12] it was 
claimed that NNSD of primes follows GOE distribution, 
while in [13, 14], the possibilities of GOE, Poisson and 
exotic Berry-Robnik [19] distribution were investigated. 
Liboff and Wong have obtained Wigner distribution and 
level repulsion for NNSD of primes by artificially includ- 
ing the gaps (no degeneracy — all primes are different) 
and 1, see [12, p. 3 113]. The gap 1 appears only once 
between 2 and 3 and should be skipped in the wake of 
infinity of primes. There is a very often reproduced figure 
showing some typical spectra (see [17, Fig. 1.2], [20, Fig. 
1.8], [21, front figure]): random levels with no correlations 
(Poisson series), sequence of prime numbers, resonance 
levels of erbium 166 nucleus, the energies a free parti- 
cle in the Sinai billiard, nontrivial zeros of the Riemann 
zeta function. In [17, p. 10] it is stated that case of prime 
numbers is far from either regularly spaced uniform series 
or the completely random Poisson series with no correla- 
tions. It is the purpose of this paper "to settle once and 
for ever that NSDD of primes follows the Poisson distri- 
bution" . The next Section II is devoted to this problem. 
In [22] M.V. Berry has calculated spectral rigidity A3 for 
zeros of the Riemann zeta function and in Sect. Ill we will 
study spectral rigidity for prime numbers. 



II. NNSD FOR PRIME NUMBERS 

In the case of primes numbers all gaps dn = Pn+i — Pn 
(except the first pair of primes pi = 2,p2 =3) are even 
integers 2, 4, 6, — These spacings are dimensionless and 
we will not perform unfolding. The unfolding obscures 
analysis of the oscillations present in the NNSD between 
original primes. Let rd{x) denote a number of pairs of 
consecutive primes smaller than a given bound x and sep- 
arated by d: 

Td{x) = ^{Pn,Pn+l < X, with Pn^l -pn = d}. (2) 



2 




FIG. 1: Plots of rd{x) for x — 2^^, 2^°, . . . , 2^^, 2^^. The histogram step widths are 2; because r2(x) ~ r4(x), therefore the visible 
step for d = 2,4 has width 4. In red exponential fits a{x)e~^^^^'^ are plotted. In the inset the plots of Td{x) / P{d) are shown. 



For odd d = 2k -\- 1 we supplement this definition by 
putting r2/c+i(x) = 0. The pairs of primes separated by 
d = 2 and d = 4 are special as they always have to be 
consecutive primes (with the exception of the pair (3,7) 
containing 5 in the middle) — in the triple of integers 
2fc + 1, 2A; + 3, 2A; + 5 the middle 2k ^3 has to be divisible 



by 3 if 2/c + 1, 2/c + 5 are prime. From the conjecture B of 
G. H. Hardy and J.E. Littlewood [23] it follows that the 
number of gaps d = 2 ("twins") is approximately equal 
to the number of gaps d = A ("cousins"), see also [6]. For 
d > 6 in [24] we have announced that 



n 

p\d,p>2 ' 



-d'K{x)/x 



for d > 6, T2(a;)( ~ T4(a;)) ^ C; 



7r2(a;) 



Co 



(3) 



Here it{x) = 0(^ — Pn) denotes the number of tegral Li(x): 
primes up to x and by the Prime Number Theorem du x x 

(PNT) is very well approximated by the logarithmic in- 7t(x) ~ Li(x) = / ~ \ ^ h . . . , (4) 

J 2 ln{u) ln{x) \n%x) 

C2 = 2np>2(l - i^^j = 1-32032... is called the 
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"twins constant". Putting in (3) 7t{x) ~ x/\n{x) the com- 
pact formula expressing rd{x) by explicitly known func- 
tions 



Td{x) ^ C2 



-d/ \n(x) 



^ ^ p\d,p>2 ^ 



(5) 



is obtained. 

During over a seven months long run of the com- 
puter program we have collected the values of Td{x) up 
to X = 2^^ ^ 2.8147 X 10^^. The data representing 
the function rd{x) were stored at values of x forming 
the geometrical progression with the ratio 2, i.e. at 
X = 2^^, 2^^, . . . , 2"^^, 2"^^. Such a choice of the interme- 
diate thresholds as powers of 2 was determined by the 
employed computer program in which the primes were 
coded as bits. The data is available for downloading 
from http: / / pracownicy.uksw.edu.pl/mwolf/gapstau.zip. 
The resulting curves are plotted in Fig.l. Characteristic 
oscillating pattern of points is caused by the product 



p{d)= n 



p\d,p>2 



p — 1 



(6) 



appearing in (3), see inset in Fig. 1. This product for the 
first time appeared in the paper of Hardy and Littlewood 
[23] and it has local maxima for d equal to the products 
of consecutive primes ( "primorials" , i.e. factorials over 
primes 2 • 3 • 5 . . . • = Pntt)- Clearly visible in Fig. 1 
are oscillations of the period 6 = 2x3 with overimposed 
higher harmonics 30 = 2 x 3 x 5 and 210 = 2x3x5x7, 
i.e. when P{d) has local maxima P(6) = 2, P(30) = 
8/3 = 2.666. . . P(210) = 16/5 = 3.2 (local minima are 
1 and they correspond to d = 2'^). We have performed 
the discrete Fourier Transform of P{d)^ i.e. we calculated 
numerically 



p(—) 

\2M) 



M-l 



,27Tkn/M 



(7) 



k=0 



where n = 0, 1, 2, . . . , M — 1 and n/2M plays the role 
of discrete frequency. Having P{f) we can calculate the 
power spectrum density /S'(/) =\P{j^)\'^ - The large value 
of S{f) at some frequency / means that the dependence 
of P{d) on d has some harmonic component of the period 
T = 1//. Thus in the Fig. 2 we have plotted S{f) versus 
l/f = dto show main periods 5, 6 = 2x3, 10 = 2x5, 14 = 
2 X 7, 30 = 2 X 3 X 5 ... of P{d). These oscillations are the 
reason why the Poisson distribution was not attributed to 
NNSD of primes in the past: e.g. P(2) = P(4) = 1 while 
P(6) = 2 and the plot should be made with logarithmic 
scale on the y axis to suppress these oscillations. 

In [25] E. Bombieri and H. Davenport have proved that: 



E n 

k=lp\k,p>2 



p — 1 



+ 0{ln\n)y, (8) 



i.e. in the limit n ^ 00 the number 2/C2 is the arithmeti- 
cal average of the product Hpi/c The main period of 
oscillations is 6 hence we can write: 



p{d)= n 



p\d,p>2 



p-1 ^ f2TTd\ 

_«« + ;3cos( — 1. 



(9) 



The numerical value of a is equal to 2/C2 to repro- 
duce the average value of P{d) in (8). It can be ex- 
plained by taking into account that cos(27r2/c/6) = 1 
while cos(27r(2A: -h 2)/6) = cos(27r(2A: -h 4)/6) -- 
hence by the equation: 



I and 



lim 



Iv-/ o /27r2fc\\ 
-^(a + /3cos^— jj=a 

k=l ^ ^ 



(10) 



the value of the parameter /3 does not contribute to 
the average of r.h.s. of (9). Thus from (8) we have 
a = 2/C2 ^ 1.5147801281. Requiring, that the com- 
bination a -h /3 cos(27r(i/6) for d = 6k takes the value 2 
times larger than for d = 6A: + 2 and d = 6k -\- 4: gives 
j3 = a/2 ^ 0.75739. Fitting of the parameters a and ^ 
can be done also numerically by standard General Lin- 
ear Least Squares, see e.g. [26]. We have used the proce- 
dure Ifit from [26] with 2500000 numbers of points: for 
d = 2, 4, 6,..., 5000000. The output of the computer run 
was: a = 1.51478, P = 0.75471 ^ I/C2. Hence we 
propose the compact formula (see inset in Fig. 2): 



p{d) 



TT PJZI 
}^ 33-2 



p\d,p>2 



(11) 



which allows to substitute for P{d) an expression more 
amenable for algebraic manipulations. Such an approx- 
imation may be relevant for calculations of correlations 
functions for zeros of the Riemann zeta function, where 
sums involving product P{d) appear very often [27]. 
Let us define the rescaled quantities: 



Td{x) = 



XTd{x) 



C2P{d)lT^{xy 



D{x,d) = 



d7T{x) 



(12) 



Because x/7r{x) ~ ln(x) is the mean distance between 
two consecutive primes up to x we see that D{x^d) corre- 
sponds to the distances between "unfolded" primes — 
normalized spacing between two consecutive primes is 
d/hi{x). From the conjecture (3) we expect that the 
points {D{x^d)^Td{x))^ d = 2,4,... should coincide 
for each x — the function Td{x) displays scaling in the 
physical terminology. In Fig. 3 we have plotted the 
points {D{x,d),Td{x)) for x = 2^8, 2^8, 2^^ If we de- 
note u = D{x^d) then all these scaled functions should 
exhibit the pure exponential decrease e~^: Poisson dis- 
tribution shown in red in Fig. 3. We have determined 
by the least square method slopes s{x) and prefactors 
a{x) of the fits a(x)e~^^^^^ to the linear parts of plots of 
{D{x,d),ln{Td{x))). The results are presented in Fig. 4. 
The slopes very slowly tend to 1: for over 6 orders of x 
they change from 1.187 to 1.136 while the prefactors drop 
from 1.512... to 1.273... . 

The explicit dependence on x of the quantities Td{x) 
and "unfolded level spacings" D{x^d) ^ d/ln(x) show 
that gaps between primes are not described by the sta- 
tionary Poisson distribution. It is a common belief that 
the Poisson NNSD of the quantum energy levels is linked 
with integrable systems with more than one degree of 
freedom. In [28] P. Crehan has shown that for any se- 
quence of energy levels obeying a certain growth law 
(l^nl < e^^+^ for some a G M+, b G M), there are 
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FIG. 2: The plot of power spectrum S{f) calculated from 
M = 2^° = 1024 values of P{d) plotted versus 1// to 
show main periods of P{d). The y axis was broken to 
make visible peaks at d 7^ 6. In the inset the plots of 
P{d) and the approximation (11) are presented. 



10^ ^ 




4 8 12 16 

D(x, d) 



FIG. 3: Plots of {D{x,d),Td{x)), [d =2,4,...) for x = 
228 ^YiQ plot of e""". Only the points 

with Td{x) > 1000 were plotted to avoid fluctuations at 
large D(x,d) due to small values of rd{x) for large d. 



infinitely many classically integrable Hamiltonians for 
which the corresponding quantum spectrum coincides 
with this sequence. Because from PNT it follows, that 
the n— th prime Pn grows like Pn ~ nln(n) the results of 
Crehan's paper can be applied and there exist classically 
integrable hamiltonians whose spectrum coincides with 
prime numbers, see also [11]. 

The smallest gap between adjacent primes is 2 (twin 
primes), while the maximal gap G{x) = maxp^^x {Pn — 
Pn-i) grows with x. We can obtain the formula for G{x) 
from (3) assuming that the largest gap up to x between 
two consecutive "levels" Pn+i — Pn appears only once: 
rG{x){x) = 1- Skipping the oscillating term P{d), which 
is very often close to 1, we get for G{x) the following 
estimation expressed directly by 7t{x): 

G(x) - (2 \n(7T(x)) - \n{x) + c) , (13) 
7r[x) 

where c = ln(C2) = 0.2778769 . . .. Substituting here the 
PNT in the form 7r(x) ^ x/\ii{x) gives the Cramer's 
conjecture [29] G{x) ~ In^(x) in the limit of large x. 
The maximal gaps G{x) are scattered chaotically, the 
largest currently known gap of 1476 follows the prime 
1425172824437699411, see [30]. The comparison of the 
above formula with real data is presented in Fig. 5. 

III. SPECTRAL RIGIDITY 

In [31] several statistical measures to describe fluctu- 
ations in the energy levels of complex systems were in- 



troduced. One which attracted much attention is the 
spectral rigidity A3. The spectral rigidity for arbitrary 
system with spectral staircase N{E) is defined as the av- 
eraged mean square deviation of the best local fit straight 
line ae-\-b to the N{E) on the interval (x, x -\- L): 



As{x; L) = ^ /min {N{x ^ e) - ae - bf de\ (14) 



The averaging procedure (•) depends on the specific 
problem, e.g. for random matrices it is the mean value 
from an ensemble of generated matrices or average over 
a set of atomic nuclei in real experiments, see e.g. [32]. 
There are in general two ways of performing the opera- 
tion mina,5, see the discussion in [31]. One can calculate 
partial derivatives of r.h.s. of (14) with respect to a and 
6, equate them to zero, solve for a, b and substitute solu- 
tions back to r.h.s. what leads to the double integrals, see 
e.g. [33, Appendix II]. These integrals can be performed 
and for the sequence £^1, £^2, • • • , of levels falling in 
the interval {x^x + L) the following explicit formula for 
As{x; L) is obtained [34]: let Ek = Ek - {x ^ L/2), then 
we have 
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FIG. 4: Plot of slopes s{x) and prefactors a(x) in 
the dependence a{x)e~^'^^^ obtained from fitting it to 
{D{d, x), ln(T,(x))) for x = 1?\ . . . , 



In the second approach the parameters a and h are ob- 
tained by fitting the straight line ax +6 to the set of points 
(xi, 7/1 ), (X2, ^2 ),..., (^n, ^n) by the least square method, 
i.e. the partial derivatives of Y^^x^Vk — CLXk — h)'^ with 
respect to a and h are calculated and put equal to zero, 
what gives the very well known expressions: 

En sr^n sr^n 

^ ^ k=l ^kVk - l^k=l l^k=l Vk 

^ELi4-(ELi^fc)' 



b = - (Vk - axk) 

In our case Xk = E^^yk = N{Ek). The spectral rigid- 
ity obtained in this second way we will distinguish from 
(15) by apostrophe A3(x; L). The formula for A'^{x; L) in 
this approach and adjusted for our problem will be given 
below, see (17). 

Spectral rigidity for primes we define by (14) with 7r(x) 
instead of N{x). To use the formula (15) the exact values 
of all primes are needed and we have used primes suffi- 
cient for calculation of /S.^{x] L) for x = 10^, 10^ and 10^^. 
As there seems to be no clear relation between the values 
of L in comparison with chosen x we have used the wide 
range of values of L: we have calculated from (15) spec- 
tral rigidity for the values of L = 2^° = 1024, . . . , 2^° = 
1.0737. . . 10^, i.e. values of L span the range over eight 



1800 




10' 10' 10' 10' 10' 10' 10' 10' 10' 10'° 10" 10'' 10'' lO'Mo'' 10'' 10" 10" 

X 

FIG. 5: The comparison of G{x) obtained from the com- 
puter search (up to x = 2^^ we have used our own data, 
for larger x we took data from the web pages [30]). For 
the plot of (13) we have used the tabulated values of 7v{x) 
available at [30]. The plot of the Cramer conjecture is 
also presented. 



(15) 

I 

orders. The results are presented in Fig. 6. It is an 
open question which values of L should be considered 
for a given x — in Fig. 6 for x = 10^ the largest L is 
larger than x while for x = 10^^ the largest L is smaller 
than X. In this approach there is no natural averaging 
procedure present in the definition (14) and in the Fig. 
6 some fluctuations are seen especially for small values 
of L. There is a crossover in the behavior of As{x; L) 
at around L* = 10^ for all three values of x and we do 
not have explanation of this behavior. For smaller L the 
dependence is roughly As{x;L) ex L^-^ while after the 
crossover there is a steeper growth of As{x] L) according 
to the relation As{x; L) cx L^-^. The points in Fig. 6 
fluctuate, especially for L below the crossover, because 
there is no averaging performed during the calculation of 
As{x; L) and number of terms in the sums in (15) is about 
n = 50 for L = 1024 and around 50000000 for L = 2^^ 
(the exact values of n are very close to the predicted by 
PNT values U{x + L) - U{x)). 

Next we will present spectral rigidity for second 
method of minimizing the r.h.s of (14) over a, 6, namely 
determination of a, b by the least square method. In the 
case of primes numbers, for large x, the smooth part of 
staircase 7t{x) given by x/\n{x) is almost linear in the 
interval {x^x-\-L) as the denominator changes from ln(x) 
to ln{x ^ L) = ln{x) + L/x -h . . . what for x > L > 1 
again is \n{x). There are a few websites [30] offering the 
tables of values of 7r{x) (as well as other number theo- 



As{x;L) 




3n 
2L2 



fc=l \fe=l / k=l 



2k + l)Ek 
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retic functions). In these data files the values of 7t{x) are 
tabulated with different step size of the best resolution 
is at the A. V. Kulsha's page: the file pi.txt of the size 
421MB contains counts of 7r(x) with a step of 10^ from 
X = 10^ to X = 2.5 X 10^^. Now we will give the formula 
for calculating the integral appearing in the definition of 
A^(x;L): 



I{L) 



f 

Jo 



{7r{x + e) 



ae ■ 



hf de 



(16) 



appropriate for our data. Let us assume, that the values 
of 7r(x) in the integral (16) are known with the resolution 
h: yk = 7r{x + kh); hence we assume that 7t{x) is constant 
on the intervals {kh, {k + l)h) (in fact 7r{x) is constant 
only between two consecutive primes). We regard this 



sampling of 7r{x) with different steps h as the averaging 
procedure hidden in the angle bracket in (14) — taking 
values of 7r{x) at all consecutive primes would introduce 
fluctuations. The combination 7r(x + e) — ae — 6 is the 
linear function on the intervals {kh, {k + l)h) and we can 
write (we assume here that L is the integer multiple of 
h): 

2 

X(L) = / {tt{x -\- e) — ae — h) de = 
Jo 



k=o -^^^ 



{yk -ae-b) de. 



Performing elementary integration we obtain: 



L/h-l 



a L 2 

A3(x;L) = 62 + a&L+^- + - ^ yk{yk - 2b)h - ayk{2k ^ l)h 



k=0 



(17) 



It should be noted, that a and b in the (17) obtained from 
fitting ae + 6 to points 7r(x + e), < e < by least-square 
method are functions of L and x, see below (19), thus we 
see that the distribution of primes is non- stationary. 

The value of A^^{x;L) given by (17) should not de- 
pend on h. To test this presumption we have calcu- 
lated A's{x;L) for xi = 10^^ and X2 = 10^^ and for 
hi = 10^/l2 = 2 X 10^/^3 = 4 X 10^. We have cho- 
sen the following sequence of values of L = 16/ii = 
1.6 X 10^^ 32hi = 3.2 x 10^^ . . . 2'^^hi = 8.388608 x 10^^ 
for both xi,X2 and additionally L = 1.5 x 10^^ for 
X2 = 10^^. It means that for h2 and hs the number of 
terms in the sum in (17) was 2^, 2^, . . . , 2^^ = 4194304 
for /i2 and 2^, 2"^, . . . , 2^1 = 2097152 for hs appropri- 
ately. For each L the a and b parameters were fitted 
by the least-square method to the points {x + kh, yk = 
7r{x + kh)), A: = 0, 1, . . . , L//i. In Figures 8 and 9 we 
present the results. Two types of behaviors are seen in 



these figures: the constant in L values of A3 depending 
on h and the collapse of plots of A3 for all h when the 
increase of A3 with L begins. It seems that to get rid 
of dependence on h the sufficiently large number L/h of 
terms in the formula (17) has to be summed up. The 
inspection of data shows, that to have the independence 
of A3 on a few thousands of terms in the sum in (17) is 
sufficient (for largest L there are millions of terms in this 
sum, see plots in royal red in Fig. 8 and 9). It is possible 
to find heuristically the values of the constant in L parts 
of A3. Namely, we consider the smooth part of 7r{x) given 
by (x + e)/ ln{x + e) and the straight line ae-\-b obtained 
by best fitting to the values of {x + kh) / \ii{x + kh). The 
experiments show, that the fits cross (x + e)/ ln(x + e) on 
the interval e G (0,L) roughly at e = L/4 and e = 3L/4, 
see Fig. 7, thus from {x + L/4)/ln(x + L/4) = aL/4 + b 
and {x + 3L/4)/ ln(x + 3L/4) = a3L/4 + 6 we get 



2 / X + 3L/4 



L14_\ 



L V Hx + 3L/4) \n{x + L/4) J ln{x) 2x In' {x) 



L 1 , , 

terms or higher 



, x + L/4 X 
b = - aL/4 



ln(x + L/4) 



ln(x) 41n'(x) 



L 1 , , 

2 h terms — or higher 



(18) 
(19) 



Using yk = {x -\- kh)/\n(x + kh) ~ (x + kh)/\n{x) — {kh)'^ /2x\ii'{x) we obtain in (17) sums over k which can be 
calculated exactly and retaining the leading terms gives: 



hL f 1 

3 + I terms — or higher 



(20) 



Because A's{x\L) > we have from above /i'/31n'(x) > /il//41n^(x), i.e. L < 4/iln(x)/3, what for x = 10^^ gi 



ives 



7 



L < 40/i. Surprisingly the first term in (20), not depend- 
ing on L but being the function of x, gives the expression 



A'six-^L; h) 



31n^(a;) 



+ .. 



(21) 



which works very weh even for L = 1024/i for xi = 10^^ 
and L = 8192/i for X2 = 10^^, as it is seen in Figures 8 
and 9, where the predicted values h? /?>\v?{x) are plotted 
by dashed lines together with the plots of ^'^{x] L; h) ob- 
tained from (17). In fact this agreement is astonishing: 
e.g. all A3 (10^^; L; hi) for initial 11 values of L have first 
three digits the same: 2.455 . . . x 10^^ while (21) predicts 
2.45588 ... X 10^^. In Fig. 8 we were able to make the plot 
for L up to almost lO^xi, while in Fig. 9 the largest L is 
smaller than ^2, thus we expect bending of A3(x2;I/; h) 
for larger L, similar to the behavior of A3(xi;L;/i) on 
Fig. 8. 

It is well known that for the stationary Poisson distri- 
bution averaged spectral rigidity is a linear function of L: 
As{L) = L/15, see e.g. [31, eq.(61)] or [33, Appendix II]. 
The case of prime numbers is non-stationary and hence 
despite the Poisson decrease P{s) = the spectral 
rigidities have more complicated behavior. Plots of both 
spectral rigidities A3 and A3 display qualitatively simi- 
lar behavior: there is a value L* above which the steeper 
increase of spectral rigidities begins and this dependence 
is L^, where 7 « 3.1 for A3 and 7 ^ 3.69 for A3. Heuris- 
tically existence of this crossover can be justified by the 
following reasoning: for moderate values of L the straight 
line ae + 6 approximates 7r{x-\-e) quite well leading to the 
small values of the integral J^^^(7r(x + e) — ae — b)^de^ 
while for larger L the discrepancy between 7r(x + e) and 
the straight line increases leading to larger values of A3's. 



IV. UNFOLDED PRIMES 

For energy spectrum E*! , , . . . one usually performs 
unfolding to pass to the dimensionless variables ei, 62, . . . 
via the definition: 



en = N{En). 



(22) 



Then the average spacing between two consecutive 
6^,6^+1 is equal to 1 and this procedure removes the 
individual properties of a system. Although primes are 
dimensionless we can perform the unfolding using the def- 
inition 



= Li(p^) 



Pn 



In(pn) " 



(23) 



Then the unfolded spacings are Dn 
ing pn^i =Pn + dn wc obtain 



Dn 



r^+i —Vn and writ- 



(24) 



and for large Pn it goes into Dn = dn/ In(pn) what agrees 
with definition of D{x^d) in (12). Because the aver- 
age distance between primes (pn-i^Pn) is In(pn) we have 



from (24) for large Pn that the average spacing between 
two consecutive (rn,rn+i) is equal to 1, as it should be. 
The values of Dn are arbitrary real numbers, while dn 
assume only even values. For example for twin primes 
Pn-\-i = Pn -\- the gap d = 2 will be mapped into 
Dn ~ 2/ ln{pn) with explicit dependence on pn. To make 
the histogram of unfolded spacings Dn the size of bin 
should be chosen. In this approach the oscillations seen 
in Fig. 1 are "smeared out" and there is no possibility 
to extract them easily from the histogram of unfolded 
gaps Dn — the behavior caused by the product P{d) is 
obscured after the change of variables d D as D de- 
pends on the value of Pn and is a continuous variable. 
To calculate spectral rigidities for unfolded "primes" 
exact values of primes Pn are needed. When the spectral 
rigidity is calculated using the values of 7r{x) at the se- 
quence of values of Xk = xo-\-kh (without knowledge of all 
primes) it is possible to define the unfolding by the map- 
ping 7r{xk) 7r{xk)\n{xk) which ensures the constant 
number of primes in intervals (x/e, x/e+i) equal roughly to 
/i, i.e. the mean level spacing will be one, as it should be. 
The plots of A3(x; L) obtained in this way are similar to 
the ones presented in Figs. 8 and 9 and we will not dwell 
on this here. 



V. CONCLUSIONS 

In this paper we have treated prime numbers as en- 
ergy levels and we applied the physical methods used to 
study spectra of quantum systems to the description of 
distribution of prime numbers. We presented large nu- 
merical data (up to X = 2.814 ... X 10^^) in support of the 
formula (3) for NNSD between consecutive primes. Ob- 
tained formulas show that the NNSD for prime numbers 
is non- stationary. It was also possible to obtain analytical 
formula (13) for the maximal difference between two adja- 
cent primes smaller than x. The case of primes numbers 
gives the rare opportunity to calculate spectral rigidity 
As{x; L) for the wide range of x and L — for real phys- 
ical systems usually only hundreds (nuclei), thousands 
or hundreds of thousands (e.g. billiards) of energies are 
known. As the main result of this paper we regard the 
scaling relations (12) and probably the first in the liter- 
ature attempt to calculate spectral rigidity As{x;L) for 
prime numbers. All the above analysis can be repeated 
for subsets of prime numbers, for example for the twin 
primes (both p and p-\-2 are prime) , cousin primes (both 
p and J9 + 4 are prime) or the primes of the form Ak^ + 1; 
in the latter case the "energy levels" are the values of k 
for which Ak^ + 1 is prime. 
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FIG. 6: Plot of A3(x;L) obtained from (15) for x = 
lO^x = 10^ and X = 10^° and L = 2^°, . . . , 2^°. Dashed 
and dotted-and-dashed lines show the character of depen- 
dencies of A3(x; L) on L before and after the cross-over. 



FIG. 7: The illustration of the experimental fact that the 
straight line best fitting (x + e)/ ln(x + e) on the interval 
e G (0, L) crosses it at e — L/A and e = 3L/4. 




FIG. 8: Plots of A^xi;L;/i) for xi = 10^^ and L - 
16/ii = 1.6 X 10^°,32/ii = 3.2 x 10^°, ... 2^2/11 = 
8.388608 X 10^^ and three values of hi = 10^, (black), 
/i2 = 2/ii (blue) and /is = 4/ii (red). On the right in regal 
red are plotted values of the number of terms L/hi — 1 
in the sum (17) and the right y axis also in regal red is 
for these numbers. The dashed lines represent values of 
(21). The coincidence of A3(xi;L)'s for all hi starts at 
approximately L = 2^^hi = 1.6384 x 10^^. 



FIG. 9: Plots of A's{x2; L;h) for xi - 10^^ and L = 
16hi = 1.6 X 10^°,32/ii = 3.2 x 10^^ , . . . 2^^ hi = 
8.388608 X 10^^ and additionally for L = 1.5 x 10^^ and 
three values of hi — 10^, (black), /i2 = 2/ii (blue) and 
/i3 = 4/11 (red). On the right in regal red are plotted 
values of the number of terms L/hi — 1 in the sum (17) 
and the right y axis also in regal red is for these num- 
bers. The dashed lines represent values of (21). The co- 
incidence of A3 (x2 ; L) 's for all hi starts at approximately 
L = 2^^hi = 2.62144 x 10^^ and follows practically power- 
like increase given by equation 2.5695 x iq-^t^^s.gss — ^^^^ 
green line presents this equation multiplied by 100, how- 
ever we expect bending of A3(x2; L; /i) for L > X2 similar 
to the one seen in Fig. 8. 



